from __future__ import division
import math
import numpy
import scipy
from scipy import integrate
from scipy.optimize import bisect
import random

##################################################################
########################## SUBFUNCTIONS ##########################
##################################################################
def isocline(x,p,n,m):
    """OVERVIEW: function that defines the isocline of perfect nestedness
    PARAMETERS:            
              x, column position in unit square
              n, the number of sites                       
              m, the number of species
              p, parameter that changes according to fill    
    RETURNS: the row postition in unit square     
    """
    y = 0.5/n+(n -1.0)/n*pow(1.0-pow(1.0-(m*x-0.5)/(m-1.0),p),1.0/p)
    return y

# ------------------------------------------------------------------------ #
def isocline_integral(p,phi,n,m):
    """OVERVIEW: calculates the residual of the integral of the isocline
                 from 0 to 1 and the fill 
    PARAMETERS: 
                 n, the number of sites 
                 m, the number of species 
                 p, parameter that changes according to fill 
                 phi, the fill
    RETURNS:     residual
    """
    parameters = (p,n,m)
    integral_out = integrate.quad(isocline,.5/m,(m-.5)/m,parameters)
    result = integral_out[0]
    return(result+phi-1)

# ------------------------------------------------------------------------ #
def get_parameter(phi,n,m):
    """OVERVIEW: gets the p value that results in the given fill
    PARAMETERS:
              phi, the fill
              n, the number of sites  
              m, the number of species
    RETURNS: p
    """
    args = (phi,n,m)
    
    try:
        parameter = bisect(isocline_integral,.001,20,args)
    except:
        if phi < 0.5: er = 'low'
        else: er = 'high'
        print 'Unable to determine the isocline because the matrix fill is to %s.' % er
        raise
        
    return parameter

# ------------------------------------------------------------------------ #
def perfnest(n,m,phi):
    """OVERVIEW: produces perfectly nested array 
    PARAMETERS:
             n, the number of sites
             m, the number of species
             phi, the fill
    RETURNS: a perfectly nested array
    """
    result = numpy.zeros( (n,m),dtype = bool )
    p = get_parameter(phi,n,m)
    for i in range(n):
        for j in range(m):
             if (n-i-.5)/n > isocline((j+.5)/m,p,n,m):
                 result[i][j] = 1
    return result

# ------------------------------------------------------------------------ #
def unexpected(n,m,p,i,j):
    """OVERVIEW: calculates the unexpectedness for a particular matrix cell
    PARAMETERS:
              n, the number of sites
              m, the number of species
              p, parameter for the isocline
              i, site index
              j, species index
    RETURNS: The unexpectedness of position i,j
    """
    unit_i = (n-i-.5)/n  #the site position in the unit square
    unit_j = (j+.5)/m #the species position in the unit square
    #calculate y and  intercept of line with slope -1 passing through unit_i,unit_j
    # if (unit_i> isocline(unit_j,p,n,m)):
        # if (cell_value==1):
            # return 0.0
    # elif (cell_value == 0):
        # return 0.0
    intercept = unit_i + unit_j
    #calculate length of line in unit square
    if intercept<=1.0:
        D =  math.sqrt(2)*intercept
    else:
        D = math.hypot(2.0-intercept, intercept - 2.0)
    #print i,j,D,unit_i,unit_j
    #arguments for optimize.bisect
    args = (n,m,p,intercept)
    #find intersection of isocline and line
    if (i,j)==(0,m-1) or (i,j) == (n-1,0):
        x_intersection = unit_j
    else:
        x_intersection = bisect(isocline_intersect,.5/m,(m-.5)/m,args)
    y_intersection = isocline(x_intersection,p,n,m)
    #calculate distance from point to isocline 
    d = math.hypot((unit_j - x_intersection),(unit_i - y_intersection))
    return ((d/D)*(d/D))

# ------------------------------------------------------------------------ #
def isocline_intersect(x,n,m,p,b):
    """OVERVIEW: a function whose zero is the intersection point, for
                 calculating unexpectedness
    PARAMETERS:
              x, column position in unit sqare
              n, number of sites
              m, number of species
              p, isocline parameter
              b, y intercept of line with slope -1
    RETURNS: the difference between the isocline and a line of slope
             -1 with y intercept of b
    """          
    return isocline(x,p,n,m) + x - b

# ------------------------------------------------------------------------ #
def unexpectedness_array(n,m,p):
    """OVERVIEW: calculates the matrix of unexpectedness
    PARAMETERS:
              n, the number of sites
              m, the number of species
              p, the isocline of perfect order parameter
    RETURNS: an nxm numpy array of unexpectedness for a given n,m, and p
    """
    unexpArr = numpy.zeros((n,m))
    for i in xrange(n):
        for j in xrange(m):
            unexpArr[i,j] = unexpected(n,m,p,i,j)
    return unexpArr

# ------------------------------------------------------------------------ #
def unexpectedness_Compare(perfect,current,i,j):
    """OVERVIEW: compares the perfectly nested array to the current array
                 to determine unexpectedness
    PARAMETERS:
              perfect, the perfectly nested array
              current, the current array
              i,j, postion to be compared
    RETURNS: True if the given position is unexpected, False otherwise
    """
    if(perfect[i][j] == current[i][j]):
      unexpectedness = False
    else:
      unexpectedness = True
    return unexpectedness

# ------------------------------------------------------------------------ #
def temperature_unpacked(inArr, phi=-1, p=-1, unexpArr=[], coldArr=[]):
    """OVERVIEW: calculates the nestedness temperature of a 2D, binary,
                 unpacked numpy array
    PARAMETERS:
              inArr, the 2D, binary numpy array to shuffle
              p, the isocline of perfect order parameter. It is calculated if not given
    RETURNS: the temperature of inArr (unpacked)
    """
    # Get some array attributes
    n,m = inArr.shape
    k = 100/0.04145 #100/Umax

    # Calculate p and phi if necessary
    if phi == -1: phi = float(inArr.sum())/(n*m)
    if p == -1: p = get_parameter(phi,n,m)
    
    # Create an array of unexpected values and a perfectly nested
    # array if necessary. This avoids having to recalculate the 
    # unexpectedness for repeat cell locations
    if coldArr == []: coldArr = perfnest(n,m,phi).astype(bool)
    if unexpArr == []: unexpArr = unexpectedness_array(n,m,p)
    inArr = inArr.astype(bool)
    
    ########## METHOD 1 ############
    # # Calculate the temperature by nested for loops (slower than below)
    # preU = 0.
    # for i in xrange(n):
        # for j in xrange(m):
            # # Check if cell is unexpected
            # if coldArr[i,j] != inArr[i,j]: preU += unexpArr[i,j]
      
    ########## METHOD 2 ############
    # Calculate the temperature using vectors (faster than above)
    I = inArr.reshape((n*m))
    C = coldArr.reshape((n*m))
    U = unexpArr.reshape((n*m))
    preU = (U[I!=C]).sum()
          
    ############ METHOD 3 ############
    ## Calculate the temperature by upper/lower matrices
    # inArr = numpy.array(inArr,dtype=bool)
    # unexpUpper = ~inArr*coldArr*unexpArr
    # unexpLower = inArr*~coldArr*unexpArr
    # U = (unexpUpper.sum() + unexpLower.sum())/(n*m)
    
    U = preU/(n*m)
    T = k*U
    
    return T
    
# ------------------------------------------------------------------------ #
def n_shuffle(vectorLen,nElems2Shuffle=-1):
    """OVERVIEW: shuffles a vector range (1D numpy array) 
    PARAMETERS:
              vectorLen, the length of the vector range to create 
              nElems2Shuffle, the number of indices to shuffle
    RETURNS: a 1xvectorLen numpy arange where nElems2Shuffle elements have been shuffled  
    NOTES:
         For nElems2Shuffle a value < 0 indicates all element indices 
         should be shuffled, and a value of 0 or 1 indicates that no element indices
         should be shuffled. If nElems2Shuffle < 0 then this returns 
         numpy.permutation(<#_elements_in_inVect>)
    """
    # Imports
    from random import sample
    
    # Get a rand ordering of element indices
    if nElems2Shuffle < 0 or nElems2Shuffle == vectorLen: shuffledVector = numpy.random.permutation(vectorLen)
    elif nElems2Shuffle in (0,1): shuffledVector = numpy.arange(vectorLen)
    elif nElems2Shuffle > 1:
        # Get the indices out to reorder and make sure they arent in order
        shuffledVector = numpy.arange(vectorLen)
        inds2Shuffle = sample(shuffledVector,nElems2Shuffle)
        while inds2Shuffle == sorted(inds2Shuffle):
            numpy.random.shuffle(inds2Shuffle)
            
        sortedInds = numpy.sort(inds2Shuffle)
        shuffledVector[sortedInds] = inds2Shuffle
    
    return shuffledVector
    
# ------------------------------------------------------------------------ #
def rearrange_array(A, newRows, newCols):
    """OVERVIEW: permutes the columns and rows of an array based on the given row
            and column permutation vectors
       PARAMETERS:
                A, the array to be rearranged
                newRows, the row permutation vector
                newCols, the column permuatation vector
       RETURNS: the rearranged numpy array
    """
    n,m = A.shape
    r1 = numpy.reshape(newRows,(n,1))
    r2 = numpy.ones(m,dtype=int)
    R = numpy.reshape(r1*r2,n*m)
    C = numpy.tile(newCols,n)
    flatIndices = R*m+C
    
    A = numpy.reshape(A,(n*m))
    return A[flatIndices].reshape((n,m))
    
# ------------------------------------------------------------------------ #
def flip_col(mat, index, n, m, p):
  """OVERVIEW: Finds the position of the isocline within the given
               column; randomly finds a cell above and a cell below
               whose values are NOT unexpected; and flips their values.
               If either side of the isocline contains only unexpected
               values, no flip is performed (because a flip would lower
               the matrix temperature).
  PARAMETERS:
              mat, a site-species presence matrix
              index, the index of the column whose cells to flip
              n, the number of sites
              m, the number of species
              p, the isocline parameter

  RETURNS:  If a flip has been performed, returns the indices of the
            affected rows.
            Otherwise, returns -1.
            (Changes to the matrix are made directly)
  """
  # To prevent flips that lower matrix temperature,
  # get positions of cells that are NOT unexpected
  cells = unflipped(mat, index, n, m, p, True)

  if(len(cells[0]) > 0 and len(cells[1]) > 0):
    toflip_above = random.sample(cells[0],1)
    toflip_below = random.sample(cells[1],1)
  else: # everything on at least one side of the isocline has been flipped
    return -1

  # transpose mat so that indices run along columns, not rows
  mat = mat.transpose()

  # flip cells
  mat[index, toflip_above] = 0
  mat[index, toflip_below] = 1

  # re-transpose to get original configuration
  mat.transpose()
  # return affected rows
  return(toflip_above, toflip_below)

# end flip_col

# ------------------------------------------------------------------------ #
def flip_row(mat, index, n, m, p):
  """OVERVIEW: Finds the position of the isocline within the given row;
               randomly finds a cell to the left and a cell to the
               right whose values are NOT unexpected; and flips their
               values.
               If either side of the isocline contains only unexpected
               values, no flip is performed (because a flip would lower
               the matrix temperature).
  PARAMETERS:
               mat, a site-species presence matrix
               index, the index of the row whose cells to flip
               n, the number of sites
               m, the number of species
               p, the isocline parameter

  RETURNS: If a flip has been performed, returns nothing.
           Otherwise, returns -1.
           (Changes to the matrix are made directly)
  """
  # get positions of cells that haven't been flipped yet
  cells = unflipped(mat, index, n, m, p, False)

  if(len(cells[0]) > 0 and len(cells[1]) > 0):
    toflip_left  = random.sample(cells[0],1)
    toflip_right = random.sample(cells[1],1)
  else: # everything on at least one side of the isocline has been flipped
    return -1

  # flip cells
  mat[index, toflip_left]  = 0
  mat[index, toflip_right] = 1

  return()

# end flip_row

# ------------------------------------------------------------------------ #
def unflipped(mat, index, n, m, p, col = False):
  """OVERVIEW: Finds all cells in the given row or column that are NOT
               unexpected.
 PARAMETERS:
            mat, a site-species presence matrix
            index, the index of the row/col whose cells to flip
            n, the number of sites
            m, the number of species
            p, the isocline parameter
            col, a boolean indicating whether to flip along columns
 RETURNS:
         A tuple containing two lists of positions of cells that
         are NOT unexpected.
         The [0] index contains cell positions above or left of the isocline.
         The [1] index contains cell positions below or right of the isocline.
  """
  # determine whether to look at row or column
  if(col == True):
    # transpose so that indices run along columns
    mat = mat.transpose()

    # get total number of cells above the isocline (i.e. to the left
    # after the transpose).
    cells_to_left = n - round(isocline((index+.5)/m,p,n,m)*n)

  else: # look across the row for total cells to left
    cells_to_left = m - round(isocline((index+.5)/n,p,m,n)*m)

  # pull out the desired row/col
  row = mat[index]

  # get unflipped (nonzero) cells to the left
  left_unflipped = numpy.nonzero(row[0:cells_to_left])[0]

  # get unflipped (zero) cells to the right
  cells_to_right = len(row) - cells_to_left
  right_unflipped = numpy.nonzero(row[cells_to_left:]-1)[0]
  # In right_unflipped, indexing starts at isocline.  Add
  # cells_to_left to indices to get the correct values
  right_unflipped = right_unflipped + cells_to_left

  # return list of indices of unflipped cells in a tuple
  return((left_unflipped, right_unflipped))

# end unflipped

# ------------------------------------------------------------------------ #
def tot_unexp(mat,n,m,p):
  """OVERVIEW: Calculates the total unexpectedness of a matrix.
  PARAMETERS:
             mat, a site-species presence matrix
             n, the number of sites
             m, the number of species
             p, the isocline parameter
  RETURNS:
         Matrix unexpectedness (a float)
  """
  # sum the unexpectedness in the cells of a matrix
  unexp = 0 # initial unexpectedness
  row_indices = set(range(m)) # the set of indices within each row
  for i in range(n):
    # Get indices of unflipped cells in row i
    unflipped_i = unflipped(mat, i, n, m, p, col = False)

    # now get indices of cells that have been flipped
    unflipped_set = set(numpy.append(unflipped_i[0], unflipped_i[1]))
    flipped_i = row_indices.difference(unflipped_set)

    # calculate the unexpectedness of flipped cells in row i
    for f in flipped_i:
      unexp = unexp + unexpected(n,m,p,i,f) # add to total unexpectedness
    # end for
  # end for
  return (1./(n*m))*unexp # scale result by matrix size
# end tot_unexp()

######################################################################    
######## MATRIX PACKING AND TEMPERATURE CALCULATION FUNCTION #########
######################################################################
def temperature(inArr,phi=-1,p=-1,coldArr=[],unexpArr=[],popSize=-1,nElite=-1,nGenerations=-1):
    """OVERVIEW: packs a 2D, binary numpy array to minimize temperature
                 using a genetic algorithm 
    PARAMETERS:
              inArr, the 2D, binary numpy array to shuffle 
              popSize, the number of solutions to create/test each generation
              nElite, the number of best solutions to keep without mutation each generation 
              nGenerations, the number of times to loop through the GA   
    RETURNS:(packed_temperature, packed_array, row_order, column_order,
             best_temp_from_each_generation)  
    NOTES: The genetic algorithm used here does not follow the one used BINMATNEST
    """
    # Define GA parameters (and other initial stuff)
    packingArr = inArr.astype(bool)
    n,m = inArr.shape
    nm = (n,m)
    nEls = n*m

    if popSize == -1: popSize = 5*int(round(math.sqrt(nEls)))
    if nElite == -1: nElite = int(0.1*popSize)
    if nGenerations == -1: nGenerations = int(750*math.log(nEls)-2500)
    
    crossRate = 0.6 # probability that two parents chosen will crossover vs create identical offspring
    mutRate = 0.1 # probability that an offspring mutates
    stopCV = 0.001 # coefficient of variation (std_dev/mean) of temperatures at which to stop the genetic algorithm
    
    # Calculate p to speed up things in the temperature calculation
    if phi == -1: phi = float(inArr.sum())/(nEls) # This will be different than the desired fill
    if p == -1: p = get_parameter(phi,n,m)
    
    # Create an array of unexpected values and a perfectly nested
    # array. This avoids having to recalculate the unexpectedness 
    # for repeat cell locations
    if coldArr == []: coldArr = perfnest(n,m,phi)
    if unexpArr == []: unexpArr = unexpectedness_array(n,m,p)
    
    # Calculate a score for rows and columns based on their sums. This
    #   is used to bias the ordering of columns and rows in proposed
    #   solutions
    cScores = packingArr.sum(axis=0).astype('uint16')
    rScores = packingArr.sum(axis=1).astype('uint16')
    
    # Create the starting population of proposed solutions
    curPop = [(numpy.arange(n),numpy.arange(m))]
    curTemps = numpy.ones(popSize)*1000.
    curTemps[0] = temperature_unpacked(packingArr,phi,p,unexpArr,coldArr)
    for i in xrange(1,popSize):
        # Create a biased-random permutation of rows and columns
        cShuffled = numpy.argsort(numpy.random.rand(m)*(cScores+1))[::-1] # [::-1] sorts descending
        rShuffled = numpy.argsort(numpy.random.rand(n)*(rScores+1))[::-1]
        
        curPop.append((rShuffled, cShuffled))
        tArr = rearrange_array(packingArr, rShuffled, cShuffled)
        arrTemp = temperature_unpacked(tArr,p,phi,unexpArr,coldArr)
        curTemps[i] = arrTemp
    
    # Sort the population ascending by temperature (descending by fitness)
    tempSort = numpy.argsort(curTemps)
    curPop = [curPop[tempSort[j]] for j in xrange(popSize)]
    curTemps = curTemps[tempSort]
        
    # Begin the genetic algorithm
    tmps = []
    for g in xrange(nGenerations):
        nextPop = []
        nextTemps = numpy.ones(popSize)*1000.
        
        # Choose the nElite best solutions
        nextPop.extend(curPop[:nElite])
        nextTemps[:nElite] = curTemps[:nElite]
        
        # Make offspring
        for i in xrange(nElite,popSize):
        
            ############ CHOOSE PARENTS ##############
            # ### METHOD 1 - Both parents from the current population ###
            # # Pick two parents to cross
            # parentInds = numpy.argsort(numpy.random.rand(popSize)*curTemps)[:2]
            # numpy.random.shuffle(parentInds) # Shuffle so the rows/columns dont always come from
                # # the better/worse of the two parents
            # parents = [curPop[parentInds[0]],curPop[parentInds[1]]]
                
            ### METHOD 2 - One parent from the current population, one "migrant" that is pseudo-randomly created ###
            # Pick one parent to cross with a randomly generated parent
            parent1Ind = numpy.argsort(numpy.random.rand(popSize)*curTemps)[0]
            parent1 = curPop[parent1Ind]
                
            # Create the other parent
            parent2 = (numpy.argsort(numpy.random.rand(n)*(rScores+1))[::-1],\
                       numpy.argsort(numpy.random.rand(m)*(cScores+1))[::-1])
            
            # Decide the ordering of parents
            parents = [parent1, parent2]
            numpy.random.shuffle(parents)
            
            ############ CROSSOVER ##############
            ### METHOD 1 - Column ordering comes from one parent, row ordering from the other ###
            # Decide if the parents cross (NOTE: Does not match BINMATNEST, but is faster)
            if numpy.random.rand() <= crossRate: # Do the cross
                offspring = [parents[0][0],parents[1][1]]
                
            else: # Offspring (before mutation) is identical to parent 1
                offspring = list(parents[0])
            
            # ### METHOD 2 - Same one as in BINMATNEST ###
            # # This matches the cross done in BINMATNEST
            # offspring = []
            # for j in (0,1):
                # p1Vector = numpy.array(curPop[parents[0]][j])
                # p2Vector = numpy.array(curPop[parents[1]][j])
                
                # offspring.append(p1Vector) # offspring is identical to parent 1 by default
                
                # # Determine if the cross happens
                # if numpy.random.rand() <= crossRate: # Do the cross

                    # # Get the possible orderings to use for the offspring
                    # n2Keep = numpy.random.random_integers(nm[j])
                    # p1Kept = p1Vector[:n2Keep]
                    # p2Uniques = numpy.lib.arraysetops.setdiff1d(p2Vector,p1Kept)
                    
                    # for s in xrange(n2Keep,nm[j]):
                        # if p2Vector[s] in offspring[j][:s]: # Take the element at random from what's left
                            # if (len(p2Uniques)) == 1: elem2Use = 0
                            # else: elem2Use = numpy.random.random_integers(len(p2Uniques)-1)
                            # offspring[j][s] = p2Uniques[elem2Use]
                            # p2Uniques = numpy.delete(p2Uniques,elem2Use)
                            
                        # else: # Take the element from parent 2
                            # offspring[j][s] = p2Vector[s]
                            # p2Uniques = numpy.delete(p2Uniques, numpy.where(p2Uniques == p2Vector[s]))
                    
            ############## MUTATION ##############
            # Mutate the offspring's row and column orderings with probability mutRate
            for j in (0,1):
                if numpy.random.rand() <= mutRate:
                    nFlips = numpy.random.random_integers(2,nm[j]) # Decides the number of elements to shuffle
                    offspring[j] = offspring[j][n_shuffle(nm[j],nFlips)]
            
            # Add the offpsring to the next generation
            nextPop.append(offspring)
            
            ############# CALCULATE FITNESS (TEMPERATURE) ##############
            # Calculate the temperature of the new offspring
            offspringArranged = rearrange_array(packingArr, offspring[0], offspring[1])
            offspringTemp = temperature_unpacked(offspringArranged,phi,p,unexpArr,coldArr)
            nextTemps[i] = offspringTemp
            
        # Update the population to the new population
        curPop = nextPop
        curTemps = nextTemps
        
        # Sort the population ascending by temperature (descending by fitness)
        tempSort = numpy.argsort(curTemps)
        curPop = [curPop[tempSort[j]] for j in xrange(popSize)]
        curTemps = curTemps[tempSort]
        
        tmps.append(curTemps[0])
        
        # End execution early if the variation in the population has dropped below a certain threshold
        if curTemps.min() == 0:
            print 'Done early on generation %i because a temperature of 0 was found.' % g
            break
        elif (curTemps.std()/curTemps.mean()) <= stopCV:
            print curTemps
            print 'Done early on generation %i because the coefficient of variation in temperatures is less than %0.4f.' % (g,stopCV)
            break
        
    # Return the array with the lowest temperature (i.e. best packed), as well as the temperature itself
    bestIndex = curTemps.argmin()
    bestTemp = curTemps[bestIndex]
    bestRows,bestColumns = curPop[bestIndex]
    bestArr = rearrange_array(packingArr, bestRows, bestColumns).astype(int)
    return (bestTemp, bestArr, bestRows, bestColumns,tmps)

#############################################################    
############## PSEUDO-DATA CREATION FUNCTIONS ###############
#############################################################
def pert_a(coldArr,unexpArr,phi,T,p):
    """OVERVIEW:  Perturbs a perfectly nested array to get an array of
               approximate desired temperature. The perturbation is based on the 
               average value of the unexpectedness array.
    PARAMETERS:
              coldArr, the perfectly nested array
              unexpArr, the unexpectedness array
              phi, the desired fill associated with the other parameters
              T, the desired output temperature
              p, the parameter that defines the isocline of perfect order  
    RETURNS: A tuple containing
             1. The pseudo-data matrix as an n by m numpy array
             2. The expected temperature of that matrix *before packing*,
                based on the the polynomial fit done in MATLAB, *if* 
                the number of cell flips wasn't then scaled by that number.
    """
    
    # Check for easy input (i.e. T = 0)
    if T == 0: return (coldArr, 0)
    
    n,m = coldArr.shape
    
    # Factor to multiply n2Shuffle by to get the proper temperature
    # The formula for this factor was generated using polyfitn (extra tool) in MATLAB
    X1 = n*m # "Size" of the matrix
    X2 = phi # fill of the matrix
    X3 = float(T) # desired output temperature
    
    expectedTOut = T # THIS NEEDS TO CHANGE!
    
    shuffleFactor = T / expectedTOut
    if shuffleFactor <= 0: shuffleFactor = 1.

    # shuffleFactor = 1;
    
    # Vectorize the inputs
    n,m = coldArr.shape
    nm = n*m
    CVec = coldArr.reshape((nm)).copy()
    UVec = unexpArr.reshape((nm)).copy()
    
    # Calculate constants
    k = 100/0.04145 #100/U_max(theoretical maximum U for any matrix
    UTarget = T/k
    uAvg = UVec.mean()
    preShuffle = shuffleFactor*((nm*UTarget)/uAvg) # number of elements to shuffle to get
        # the desired temperature
    if preShuffle > nm: preShuffle = nm
    
    # Get the indices of zeros and ones to flip
    n2Shuffle = int(round(preShuffle/2.))
    trueInds = numpy.random.permutation(numpy.nonzero(CVec==True)[0])[:n2Shuffle]
    falseInds = numpy.random.permutation(numpy.nonzero(CVec==False)[0])[:n2Shuffle]
        
    # Do the perturbation
    COrder = numpy.arange(nm)
    COrder[trueInds] = falseInds
    COrder[falseInds] = trueInds
    CVec[trueInds] = False
    CVec[falseInds] = True
    return (CVec.reshape((n,m)), expectedTOut)

# ------------------------------------------------------------------------ #
def pert_n(coldMat,unexpMat,phi,T,p):
    """OVERVIEW: produces a random array according to Nick's Method
    PARAMETERS:
              n, the number of sites
              m, the number of species
              phi, the fill
              T, the temperature
    RETURNS: a tuple containing
             1. A randomly nested array via Nick's Method 
             2. The expected temperature of the array *before* packing
    """
    k = 100/0.04145 #100/U_max(theoretical maximum U for any matrix
    U_target = T/k
    U_current = 0.0
    # p = get_parameter(phi,n,m)
    n,m = coldMat.shape
    while U_current<U_target:
        i = numpy.random.randint(0,n)
        j = numpy.random.randint(0,m)
        # u = unexpected(n,m,p,i,j)
        u = unexpMat[i,j]
        if u == 0:
            if coldMat[i][j]==1:
                coldMat[i][j] = 0
            else:
                coldMat[i][j] = 1
            # u = unexpected(n,m,p,i,j)
            u = unexpMat[i,j]
            U_current = U_current + u/(n*m)
    return (result, U_current)

# ------------------------------------------------------------------------ #
def pert_c(n,m,phi,T,mat='',p='', err = 1, max_iter = 100):
  """OVERVIEW: Creates a site-species presence matrix with a given fill
               and temperature.
  PARAMETERS:
             n, the number of sites
             m, the number of species
             phi, the matrix fill, i.e. the proportion of nonzero elements
             T, the desired temperature of the matrix
             err, half-width of acceptable interval around T
                  e.g. err = 1 (default) stops function execution
                  when matrix temperature is within 1 degree of T
             max_iter, maximum iterations before function terminates
  RETURNS:
         A tuple containing the matrix and its actual temperature.
  NOTES:
         Some combinations of T and phi cause the algorithm
         to converge on a lower temperature than desired.
  """
  # First, create a perfectly nested matrix with dimensions n,m and fill phi
  if mat == '': mat = perfnest(n,m,phi)
  # keep track of successive versions of mat.  If the algorithm overshoots temperature T,
  # it will step back to a previous version of the matrix and try again.
  mat_versions = [mat.copy()]
  
  # Do the easy check:
  if T == 0: return (mat, 0)

  if p == '': p = get_parameter(phi,n,m) # get isocline parameter
  temp = 0  # current temperature of mat
  prev_temps = [temp] # keep track of previous temperatures, to check for convergence
  # get unexpectedness associated with temperature T
  # T = k*U, where k = 100/Umax = 100/0.04145
  target_unexp = T / (100/0.04145)

  # counter and stopping condition
  i = 0
  done = False

  # Iteratively perturb and pack the matrix until either T or max_iter is reached
  while done == False and i < max_iter:  # loop level 1
    i = i + 1
    # print i
    # calculate current unexpectedness
    unexp = tot_unexp(mat, n, m, p)
    cols_done = False # columns contain cells to be flipped

    # perturb matrix by flipping cells until either target_unexp
    # is reached or no more cells are left to flip.
    j = 0
    while unexp < target_unexp and cols_done == False: # loop level 2
      j = j + 1

      col_indices = range(m) # all possible columns
      need_new_index = True
      # Randomly select columns to perturb until either a successful flip occurs
      # or no more columns are left to try
      while need_new_index == True:  # loop level 3
        col_index = random.sample(col_indices,1)[0]
        col_flipped = flip_col(mat, col_index, n, m, p)

        if col_flipped != -1: # flip was successful
          need_new_index = False
        else: # no flip; one side of isocline has been completely flipped already
          if len(col_indices) > 1:
            col_indices.remove(col_index) # remove col_index from sample pool
          else:
            col_indices.remove(col_index)
            need_new_index = False        # no more columns left to choose
      # end loop level 3

      # If all columns have been tried unsuccessfully, break out of loop level 2,
      # i.e. go directly to matrix packing and temperature calculation.
      if len(col_indices) == 0:
        cols_done = True
        continue
      else: # flip cells in rows that were affected by the column flip; get unexpectedness
        row_flipped = flip_row(mat, col_flipped[0][0], n, m, p)
        row_flipped = flip_row(mat, col_flipped[1][0], n, m, p)

        # get total unexpectedness = sum of cells in the unexpectedness matrix
        unexp = tot_unexp(mat,n,m,p)

        if unexp > target_unexp:  # more than enough unexpectedness. break out of loop level 2
          cols_done = True
          continue
        # else: keep perturbing the matrix
    # end loop level 2

    # Now the matrix has either reached target_unexp or cannot be further perturbed.
    # pack the matrix and get its temperature.
    temp_info = temperature(mat)
    temp = temp_info[0]
    mat  = temp_info[1]

    # Check for early convergence and stopping conditions
    if temp == prev_temps[-1]: # matrix has converged on a temperature less than T
      # step back to the previous matrix
      mat = mat_versions[-1].copy()
      # discard most recent matrix unless only the first matrix remains
      if len(mat_versions) > 1:
        mat_versions = mat_versions[:-1]
      #prev_temps   = prev_temps[:-1]
      # don't discard most recent temperature- otherwise the algorithm won't recognize
      # early convergence
    elif temp < T-err: # matrix temperature is below T
      # update mat_versions and prev_temps
      mat_versions.append(mat.copy())
      prev_temps.append(temp)
    elif temp > (T-err) and temp < (T+err):  # matrix temp is close enough to T. terminate loop.
      done = True
    else: # temp > T+err; matrix temperature is too high.
      # step back to previous matrix (unless the current matrix is the first matrix)
      if len(mat_versions) > 1 and len(prev_temps) > 1:
        mat = mat_versions[-1].copy()
        mat_versions = mat_versions[:-1]
        prev_temps   = prev_temps[:-1]

    if i == max_iter:  # user-defined maximum iterations have been reached
        # return the matrix with the highest unexpectedness
        # print "choosing best matrix..."
        unexps = []
        if len(mat_versions) > 1:
          for mv in mat_versions:
            unexps.append(tot_unexp(mv,n,m,p))
          best_unexp = max(unexps)
          best_index = unexps.index(best_unexp)
          mat = mat_versions[best_index]
          temp =  temperature(mat)[0]
        # else, the only previous matrix to choose from is the perfectly nested one.
        # stick with the most recent matrix.

  # end loop level 1

  return (mat, temp)
# end pert_c()

# ------------------------------------------------------------------------ #
def generate(n,m,phi,T,method='pert_a',calcTemp=True):
    """OVERVIEW: Uses one of several methods included in this file to   
               create pseudo-data. Also calculates the output matrix's
               temperature if desired.
    PARAMETERS:
              n,m, the number of rows and columns in the matrix
              phi, the desired fill of the output matrix    
              T, the desired temperature of the output matrix (a seed value) 
              method, the method used to create the pseudo data. Can be any 
                      in {pert_a, pert_n, pert_c}. Default is pert_a
              calcTemp, a boolean to indicate whether or not the actual output of 
                        matrix temperature should be calculated. default is true
    RETURNS: a tuple (outputArr, TOutExpected, TOut)
    """             
    # Define constant inputs
    coldArr = perfnest(n,m,phi)
    p = get_parameter(phi,n,m)
    unexpArr = unexpectedness_array(n,m,p)   
    
    # Create the pseudo data using available methods
    if method == 'pert_a':
        outputArr, TOutExpected = pert_a(coldArr,unexpArr,phi,T,p)
        
    elif method == 'pert_n':
        outputArr, TOutExpected = pert_n(coldArr,unexpArr,phi,T,p)
        
    elif method == 'pert_c': # Placeholder for Christine's code
        outputArr, TOut = pert_c(n,m,phi,T,coldArr,max_iter=10)
        TOutExpected = T
        
    # Calculate the temperature of the output matrix, if indicated
    if method == 'pert_c': pass
    elif calcTemp: TOut,outputArr = temperature(outputArr,phi=phi,p=p,coldArr=coldArr,unexpArr=unexpArr)[0:2]
    else: TOut = -1
        
    return (outputArr.astype(int), TOutExpected, TOut)